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Abstract 

We have measured the specific heat of zincblende ZnS for several isotopic compositions and 
over a broad temperature range (3 to 1100 K). We have compared these results with calculations 
based on ab initio electronic band structures, performed using both LDA and GGA exchange- 
correlation functionals. We have compared the lattice dynamics obtained in this manner with 
experimental data and have calculated the one-phonon and two-phonon densities of states. We 
have also calculated mode Griineisen parameters at a number of high symmetry points of the 
Brillouin zone. The electronic part of our calculations has been used to investigate the effect of the 
3d core electrons of zinc on the spin-orbit splitting of the top valence bands. The effect of these 
core electrons on the band structure of the rock salt modification of ZnS is also discussed. 

PACS numbers: 63.20.-e, 63.20.dk, 63.20.D-, 68.35.bg, 65.40.Ba, 71.55.Gs, 71.70.Ej 
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I. INTRODUCTION 



This article is part of an on-going effort to investigate, both theoretically and exper- 
imentally, the specific heat of semiconductors and semimetals, with emphasis in the low 
temperature region where strong deviations from the Debye law take place.-i^i^ The 



reader will find more references to our recent work on the subject in Refs. |^ andlsj. Most of 
our past experimental and theoretical work has been performed for materials with different 
isotopic compositions, the exceptions being the mercury chalcogenides (HgX, X=S, Se, Te) 
for which only natural elements were considered. For the HgX materials, however, we made 
a theoretical investigation of the effect of spin-orbit (s-o) interaction on the specific heat and 
the phonon dispersion relations (this investigation was also carried out for Sb,-, Bi,-, and 
the lead chalcogenides^). The ab initio calculations of Ci,, based on the electronic structure 
as described by local density functionals, are able to account rather well for the experi- 
mental results.- The availability of relativistic electronic band structures as a by-product of 
our calculations, suggested their use to investigate a number of interesting electronic and 
structural properties such as s-o splittings (including the corresponding linear terms of the 
valence bands around k=0) , and also phase transitions under pressure.- 

In this article we extend this work to ZnS, a material for which several isotopes of its 
constituent atoms and single crystals of its zincblende variety were available to us.- 

ZnS is found in nature mostly in the zincblende (also called sphalerite, 3C) cubic structure 
(it actually gives its name to this structure-type) and constitutes the standard ore for the 
smelting of zinc. ZnS is also found as the mineral wurtzite, which lends its name to the 
canonical 2H structure in which several tetrahedral semiconductors crystallize (e. g. GaN, 
CdS, CdSe). As we shall see in the course of this work, the enthalpies of these two crystalline 
varieties of ZnS differ only very little, the 3C modification being slightly more stable than 
wurtzite except at temperatures higher than ~1200 K.— A transition from the 3C (zb) 
phase to the rocksalt (rs) structure takes place under pressure at about 15 GPa.^^ The zb 
and wurtzite (wz) phases of ZnS differ only in the stacking order of crystal planes along one of 
the three-fold symmetry axes, a fact which probably accounts for the small difference in their 
enthalpies. A large number of similar phases (more than 200), with related stacking orders 
but larger primitive cells have been identified.— They provide one of the most conspicuous 
examples of the phenomenon called polytypisnt^ We shall not be concerned in this article 
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with this fascinating aspect of ZnS. From the practical point of view, we would also like 
to mention that polycrystalline ZnS is available commercially as windows and domes for 
infrared radiation detection devices.— Because of its large electronic band gap, it is also 
useful for optoelectronic applications involving the short wavelength part of the optical 
spectrunt^ii^ although it has recently found a strong competitor in GaN.— 

We present in this paper measurements and ab initio calculations of the specific heat 
Cp of natural and isotopically modified zb-ZnS ( ®^Zn, ^^Zn, '^^S, ^^S) between 3 and ~200 
K (for ZnS with the natural isotope composition the work has been extended to 1100 K). 
We pay particular attention to the maximum in Cp/T^ which occurs at about 22 K and 
obtain from the experimental data the logarithmic derivatives of Cp vs. the isotopic masses 
of Zn and S. These results are compared with the temperature dependence of Cp in the 
manner followed in our earlier publications.- Since the calculation of Cy requires that of 
the phonon dispersion relations, we compare the results obtained for the latter from ab 
initio electronic band structures with experimental results obtained by inelastic neutron 
scattering. We also present the calculated one-phonon density of states (DOS), explicitly 
separating the vibrational components of the two constituent atoms and provide the densities 
of two-phonon pairs with a total zero wavevector, including both sums and differences of 
two phonons. These DOS are relevant to the interpretation of Raman and ir spectra. We 
compare the results obtained with extant second order Raman scattering data. 

We have performed the ab initio electronic calculations with and without s-o interaction. 
Because of the smallness of the s-o effects, they barely influence the vibrational properties, 
contrary to what was found in materials containing heavier atoms.-'^'^'^ It was found in Ref. 
Islthat the s-o splitting of core d electrons of the cation provides a negative contribution to the 
s-o splitting of the top valence bands (at the F point) which, in the case of HgS, even reverses 
the sign of this splitting.— We have investigated the possibility of such an effect in zb-ZnS and 
found that, while the s-o splitting remains positive, the admixture of 3d electrons of the Zn- 
core reduces its value. This reduction depends strongly on pressure. In the rs modification, 
however, parity prevents the 3(i-zinc admixture with the 3p wavefunctions of sulfur and 
such effects do not occur. In order to ascertain the consequences of this fact, calculations 
for rs-ZnS were performed at a pressure slightly above that of the phase transition. We 
noticed that the lack of p-d admixture at the top of the valence band is also responsible for 
converting the rs phase of ZnS into an indirect gap semiconductor (the zb gap is direct) in 
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a way similar to that found earlier for AgBr and AgCl.— The lack of inversion symmetry in 
zb also generates s-o splittings linear in k at the top of the valence bands. We evaluated 
them and compared them with previous semiempirical results. 

Since some of the work just described requires calculations of the electronic structure 
as a function of pressure, we also evaluated the phonon dispersion relations of zb-ZnS at 
several different pressures below 15 GPa. We calculated the mode Griineisen parameters 
at the r, X and L point of the Brillouin zone. From the parameters of the LO and TO 
modes at F we evaluated the Griineisen parameter of the dynamical charge and found it to 
be negative, like for most other tetrahedral semiconductors (exception: SiG^), a fact that 
is usually interpreted as signaling the decrease of ionicity with increasing pressure. 

This article is structured as follows: In Section [III we discuss the numerical procedures 
and codes employed for the ab initio calculations. Section IIIII presents and discusses some 
details of the electronic band structure calculations that were used in the subsequent sections. 
Particular attention is paid to the effects of s-o interaction mentioned above and also to the 
optimized lattice constants and bulk moduli which are obtained with the various codes 
employed. These results are compared with experimental data. 

Section [IV] is devoted to the ab initio calculations of the phonon dispersion relations 
and comparison with INS and Raman data. For the purpose of calculating Cy{T) we need 
the one-phonon DOS which we present here together with its decomposition into Zn and 
S contributions. We also use it to estimate the contribution of isotopic disorder to the 
Raman linewidths for natural and for ZnS composed of (^^Zuq.s ^^Zuq.s) (^"Sq.s ^"^Sq.s) isotopic 
mixtures. We also present the calculated two-phonon optical (k=0) DOS, for sums and for 
difference modes, and compare it with extant Raman data. 

Although some information is already available in the literature, we evaluate in Sec JVIl 
the T=0 enthalpies of the three most important phases of ZnS vs. pressure. The results 
are presented as differences with respect to that of the most stable zb phase. At T = 
the enthalpy difference of the wz phase is basically independent of pressure (in recent work, 
this difference has been shown to decrease with increasing T and to vanish for T ~ 1200 
K— ). The zb to rs phase transition is calculated to take place for T ~ between 15.8 
and 16.7 GPa, depending on whether the LDA or the GGA approximation is used for the 
exchange and correlation potential. As a by-product we report the pressure dependence 
of the phonon modes at F, L and X, and the corresponding Griineisen parameters. From 
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those parameters for the LO and TO phonons at F, we derive the volume dependence of the 
dynamical effective charge and the corresponding Griineisen parameter. Calculations of the 
electronic band structure vs. pressure reveal an anomalous decrease of the s-o interaction 
with increasing pressure, related to an increasing admixture of Zn-3(i core wavefunctions 
to those predominantly S-3p-like at the top of the valence bands. Section IIVI presents the 
measurements of Cp{T) for natural ZnS and for several isotopically modified samples. The 
results for natural ZnS are compared with various sets of rather incomplete data available 
in the literature and with our ab initio calculations. The logarithmic derivatives of Cp 
vs. isotopic mass are also compared with the derivative of Cp vs. T following the scheme 
used earlier for other binary semiconductors.-!^ Section IVIII presents the conclusions and 
summarizes the results. 



II. THEORETICAL DETAILS 

Semi-empirical calculations of the band structure and the lattice dynamics of 3C-ZnS have 
been available for quite some time.-i^ii^^ In the past few years, a number of articles reporting 
ab initio studies of electronic, structural, and vibronic properties have appeared.— i^i^ These 



articles cover some partial aspects of the work described here. Reference [id, for instance, 
reports calculated values of Cy only at five temperatures (300, 600, 900, 1200, and 1500 K) 
and thus misses the interesting region from 10 to 50 K in which the maximum of C^/T^ 
takes place. We thus decided to perform independent similar ab initio calculations of all 
aspects related to the present work. The results will be compared with those of the Chinese 
researchers when appropriate. 

The calculations were carried out using two different implementations of density func- 
tional theoryS^i^^, in order to ensure that our conclusions are independent of the basis wave 
expansion. We have also used different exchange-correlation functionals (LDA,GGA) for the 
same reason. 

We have utilized a linear response approacb2Ii2^i2£ together with an iterative minimiza- 
tion norm-conserving pseudopotential plane-wave method as implemented in the ABINIT 
package.— The local-density approximation (LDA) was employed for the exchange and cor- 
relation energy.— These pseudopotentials are single projector, ordinary norm conserving, 
based on the TrouUier- Mart ins method.— "^i^ 12 and 6 "valence" electrons were used for Zn 
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and S, respectively. A 40 Ry cutoff was defined for the plane wave expansion and a 6x6x6 
regular shifted k-point mesh. The calculation is well converged: increasing the plane-wave 
cutoff to 50 Ry and the k-points mesh to 8x8x8 yields, on average, a change of less then 
two percent in the phonon frequencies. For the phonon frequencies, we have used a q-mesh 
of 12x12x12, to guarantee a good coverage of the dispersion relations. The dynamical 
matrices were obtained from perturbation theory^i2£ and a Fourier interpolation was used 
to increase the mesh sampling to improve the description of quantities such as vibrational 
density of states, heat capacity, etc.. 

In addition we have also performed calculations with the Vienna ab initio simulation 



by; 



using the pseu- 



35 and references 



package (VASP), where density functional theory has been implemented 
dopotential method to describe the electron wave functions (see Ref. 
therein). The exchange-correlation energy was either taken in the local density approxi- 
mation (LDA) with the Ceperley-Alder prescription^!^ or with the generalized gradient 
approximation, GGA, with the PBE prescription.— The projector- augmented wave (PAW) 
scheme^ was adopted and the Zn 3d orbitals were included explicitly in the calculations. The 
set of plane waves extended up to a kinetic energy cutoff of 380 eV. This cutoff was found 
necessary to achieve highly converged results within the PAW scheme. We used a dense 
Monkhorst-Pack grid for the Brillouin zone (BZ) integration to ensure highly converged re- 
sults (to about 1-2 meV per formula unit). We also employed an accurate algorithm during 
the calculations in order to assure very well converged forces in the calculation of the dynam- 
ical matrix. At each selected volume, the structures were fully relaxed to their equilibrium 
configuration through the calculation of the forces on atoms and the stress tensor.— In the 
relaxed equilibrium configuration, the forces are less than 0.002 eV/A and the deviation of 
the stress tensor from a diagonal hydrostatic form is less than 0.1 GPa. 

In the case of VASP, highly converged forces are required for the calculation of the 
dynamical matrix using the direct force constant approach (or supercell method).— The 
construction of the dynamical matrix at the F point of the BZ is particularly simple and 
involves separate calculations of the forces in which a fixed displacement from the equi- 
librium configuration of the atoms within the primitive unit cell is considered. Symmetry 
arguments were taken into account for reducing the number of such independent distortions. 
Diagonalization of the dynamical matrix provides both the frequencies of the normal modes 
and their polarization vectors. It allowed us to identify the irreducible representation and 
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the character of the phonons. Supercell calculations were carried out in order to obtain 
phonon dispersion curves. We start with a 2x2x2 supercell for which the direct method 
delivers correct phonon frequencies at special points of the BZ which are compatible within 
the selected supercell size. If the interaction range would cease within the selected supercell 
size the phonon frequencies would provide good interpolation values for all the wave vec- 
tors. (The LO/TO splitting cannot be included using the direct method, but the dynamical 
matrix can be supplemented with a non-analytical term, which depends on the Born ef- 
fective charge tensors and the electronic dielectric constant, in order to obtain the LO/TO 
splitting.— Thermodynamic properties of the crystals at constant volume are determined by 
phonons.— The ones measured at constant pressure require an anharmonic correction which 
has been included by hand.- 

We display in Table [T] the lattice parameters we have used in our calculations which were 
found through energy minimization. The results presented below were obtained with either 
the ABINIT code (local density functional, LDA) or the VASP code (using either the LDA 
functional, or the generalized gradient approximation (GGA)). Note that the experimental 
data for ao lie between the GGA and the LDA results, the former being higher than the 
measured ones, the latter lower, as is usually the case.— 

TABLE L Lattice parameters ao and first (Bq), second (Bq) and third order {B") bulk moduli at p 
= T = 0, obtained by energy minimization with the VASP code using GGA local density functional 
( with and without s-o interaction) and the LDA (with s-o) interaction. Calculated values were 
obtained through fits to the Birch-Murnaghan equation.— The experimental values are taken from 
Ref ' 
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. results of our ABINIT-LDA values are also listed. 





ao (A) 


Bo (GPa) 


B'o 


B'^ (GPa-i) 


VASP-GGA, NOSO 


5.450 


71.12 


3.88 


0.29x10"'' 


VASP-GGA, SO 


5.450 


70.51 


3.98 


-0.85xl0~3 


VASP-LDA, SO 


5.302 


90.12 


3.79 


-0.17xl0~3 


ABINIT-LDA, NOSO 


5.319 


85.86 


4.48 


-0.86x10-3 


experiment 


5.404 


76.9 


4.4 





Correspondingly, the calculated bulk modulus Bq is larger for the LDA than for the 
GGA calculations. The higher order bulk moduli B' and B" were obtained by fitting total 



energy calculations vs. the lattice parameter qq to the Birch- Murnaghan equation.— B' is 
found to have the typical values close to 4. B" is rather small and can be neglected for all 
practical purposes. Because of our previous interest in s-o coupling effects on thermodynamic 
properties^ we have actually performed GGA calculations with and without s-o-interaction. 
This interaction does not alter the value of although it slightly increases i?o, i-e., the 
average acoustic phonon frequencies. This effect, however, is too small to be of interest. 



III. ELECTRONIC BAND STRUCTURE 

Figure [1] displays the electronic band structure of zb-ZnS as calculated with the VASP 
code using the LDA exchange-correlation Hamiltonian. Figure [T](a) shows the region around 
the direct energy gap at F as well as the bands which correspond to the "id " semi-core" elec- 
trons of Zn (between 6 and 7 eV below the top of the valence band), which were explicitly 
included in the calculation. These bands are responsible for some of the interesting phe- 
nomena to be discussed below (the small s-o splitting at the top of the valence bands and 
its anomalous dependence on pressure, the conversion of the direct into an indirect gap in 
the rs phase). 

The calculated direct gap Eq (2.19 eV, see Fig. [l](b)) is considerably smaller than the one 



measured at 19K (3.85 eV, cf. Ref. l44l ). This fact applies to all LDA and GGA calculations 
and is usually referred to as the "gap problem". It has been attributed to a discontinuity in 
the local density functional which takes place when crossing the gap.— Recent calculations 
using the GW approximation yield for this discontinuity an energy of 2.05 eV.— When 
adding this discontinuity to the 2.2 value of the gap (Fig. [U^b)) one obtains 4.25 eV. This 
energy should be lowered by the electron phonon interaction about 0.1 eV.— The resulting 
value of 4.25 eV is somewhat higher than the experimental one. Smaller values of the 
discontinuity have been very recently suggested.— The use of the experimental oq instead of 
the one obtained through energy optimization would also lower the calculated gap and thus 
decrease the remaining discrepancy. We have also performed calculations using the GGA 
for the exchange- correlation functional, also including s-o interaction. The gap obtained is 
2.07 eV, somewhat smaller than the LDA gap. The difference (0.12 eV) is less than that 
estimated from the difference in lattice constants for the two functionals (see Table I) using 
the measured value for the pressure coefficient of Eq [dEo/dp = 0.057 eV/GPa, estimated 
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FIG. 1: (a) Electronic band structure of ZnS calculated with the VASP code using an LDA exchange 
correlation Hamiltonian. (a) includes the 3d core-like electron bands of Zn. (b) Magnified part of 
the band structure of ZnS which concentrates on the bands around the energy gap, so as to display 
the s-o splittings. 

LDA-GGA gap difference 0.37 eV). 

We discuss next the value of the s-o splitting Aq of the top valence bands at F. It was 



found in the calcu 
= 0.067 eV, Ref. 



ations of Fig. [T]to be 0.062 eV, in agreement with the measured value (Aq 
50l ) and with the ab initio calculations of Carrier and Wei (0.066 eV).— 



This value is rather small compared with estimates in the literature based on the splittings 
of the 3p atomic levels of Zn and those of S (0.097 eV, Ref. |52| ). This fact suggests a 
significant contribution of the 3d core electrons of Zn, similar to that which accounts for the 
.g„ .eve.sal of A. ,„ HgS. CuC. Z„0 and CuGaS, (Ref. In oMe. to venfy this con.ect.e, 
and to obtain a quantitative estimate of the Zn-Srf contribution to Aq , we shall present in 
Sec. |Vl] calculations of the electronic band structure of rs-ZnS. The crystal structure of 
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FIG. 2: Splitting of the Fg valence bands for k along [100] and [111]. The fitted linear coefficient 
Ck is given next to the corresponding curve. Their average is -20.7 meVxA. Note that for k along 
[111] two of the quadruplet bands do not split. The splitting is itCfc • /c for k along [100] and 
±V2Ck • /c for k along [111]. 

rocksalt is centrosymmetric and thus no p-d hybridization is possible at the F point of the 
Brillouin zone. Consequently no compensation of the 3p splitting at F by the 3d splitting of 
the Zn core levels should take place. We calculate for the rs phase a splitting of 0.112 eV, 



which corresponds to that listed in Ref. 



52 



Another interesting effect of the spin-orbit interaction is the appearance of linear terms 
in k in the valence bands around F.-i^i^ The existence of these terms follows from the lack 
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of inversion symmetry through the couphng of the Fg valence band states with the Zn-Sd 



541 and a semiempirical 



"semi-core" states via the k-p and the s-o Hamiltonians. In Refs. 
perturbation expression for the coefficient is given. For the zb II- VI compounds it reduces 
to: 



Cfc 



-350- 



(in meV x A) 



(1) 



Arf c is the s-o splitting of the outmost cation d levels and E^^c their energies, respectively. 

Using for the energy denominator of Eq. ([1]) the value of 6.5 eV extracted from Fig. [1] 
and the atomic Zn-3d splitting of 0.4 eV (Table X in Ref. [sj) we find Ck = -21.5 meVxA.-^ 
This value of Ck is slightly larger than that listed in Table X of Ref. l5J for ZnS {Ck = - 15.5 
meVxA) because a smaller value (6.5 eV) of the denominator of Eq. ([T]) has been used (9 



eV was used in Ref. 



Our band structure calculations also enable us to obtain ab initio the (linear) splittings 
under consideration. The corresponding results are shown in Figure [2] for k along [100] 
and [111]. Along [100] the Fg bands split into two doublets. Along [111] Fg splits into two 
singlets, and a doublet which is not affected by the s-o interaction. From the initial slopes 
of the curves in Fig. [2](b), an average Ck = -20.7 meVxA is obtained. 



IV. LATTICE DYNAMICS 

Figure [3] displays the phonon dispersion relations of rs-ZnS, calculated using the ABINIT 
code with the LDA approximation. We have also plotted in this figure the phonon fre- 
quencies measured at room temperature by Vagelatos et al. with INS.— Because of the 
anharmonic temperature shift the measured frequencies are expected to be somewhat lower 
than calculated ones, a fact which seems to be supported by Fig. O— This difference, how- 
ever, may be enhanced by the fact that the LDA value of ao is smaller than the experimental 
one (Table [T]). We have also included in Fig. [3] frequencies calculated with the CASTEP code 
using the GGA.— These frequencies are also lower than those obtained with the LDA, even 
lower than those obtained experimentally, a fact that suggests that the difference is mainly 
due to differences in the lattice parameters (larger for the GGA than either the experimental 
ones or those calculated with the LDA, Table [11). Note that the dispersion relations of the 
LO phonons are almost flat whereas those of the TO counterparts show a slight upward 
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bending starting at the T point. This bending, however, is much less than that found for 
HgS (Ref. [s]), as corresponds to the smaller difference between the cation and anion masses. 
The data of Fig. [3] at the F point allow us to calculate the transverse (or Born) effective 
charge e* using the standard expression: 

2 2 _ 4vr (e* eY 

where wlo and cjto are the corresponding phonon frequencies at the F point (Raman 
frequencies), M the reduced mass of the primitive cell (21.52 amu for natural ZnS), Vc the 
volume of that cell, eoo the electronic contribution to the long wavelength dielectric constant, 
and e the elementary charge. Using the LDA values of the Vc (flo/4, as obtained from Table 



[T]), = 5.2 (Ref. 14^, and the calculated Raman frequencies we obtain e* = 1.76, which 
agrees reasonably with the values obtained from the experimental frequencies (2.01). In Sec. 
IVll we shall discuss the dependence of this charge on pressure. 



400 




L r X r 



FIG. 3: (color online) Phonon dispersion relations of zb-ZnS, calculated with the ABINIT-LDA 

561 ). The (green) squares represent the Raman 



code, compared with the results of INS (dots, Ref. 
frequencies measured at low temperatures. The stars display ab initio calculations performed with 
the CASTEP-GGA code^ 

In order to calculate the temperature dependence of the specific heat we need the phonon 
DOS, which can be obtained by standard procedures from the dispersion relations of Fig. 
131 We display in Fig. H] the total phonon DOS calculated for the ABINIT-LDA dispersion 
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relations, together with the separate contributions of the Zn and S vibrations. As expected, 
the acoustic branches are mainly Zn-like whereas the optic ones are mainly S-like. These 
separate components are of interest for an estimate of the broadening induced by isotopic 
fluctuations in natural crystals or in those grown with strong isotopic disorder .-i^ii^ For this 
broadening to occur, there has to be significant phonon DOS at the frequency of interest, 
projected on the atom whose disorder is being considered. This effect is shown to be impor- 
tant for ZnSe with regards to the isotopic disorder of Zn.— In Ref. o] an unsuccessful attempt 
was made to identify such effect for the Raman phonons of ZnS. This can be understood on 
the basis of the disperson relations of Fig. [3] and the DOS of Fig. HI At the frequency of the 
TO phonons at F the DOS vanishes. This is not so easy to conclude for the nearly fiat LO 
branch. However, the experimental results of Serrano et al. (their Table III) suggest that 
isotopic disorder is absent for the LO (7) phonons of ZnS.- 
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FIG. 4: (color online) Phonon DOS of ZnS obtained from the dispersion relations shown in Fig. 
[3l The DOS has been normalized so as to correspond to a total of 6 states per primitive cell. The 
partial components of the DOS (PDOS) corresponding to the Zn and the S vibrations are explicitly 
shown (red) dotted line: Zn PDOS; (blue) dashed line: S PDOS). Notice the large gap between 
acoustic and optic phonons. 

Phonon DOS have been obtained previously with semiempirical techniques.— 1^ They 
agree semiquantitatively with the ab initio results reported here and with the recent ab 
initio work of Yu et al..— These authors, however, did not report the partial components of 
the DOS. 
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While several calculations of the one-phonon DOS are already available, the spectral 
distribution of two-phonon (2Ph) DOS (with ki + k2 = 0, as imposed by optical absorption 
and Raman scattering) is not that readily available. Semi-empirical calculations of the 2- 
phonon DOS are presented in Ref. |2J in the form of histograms which have poor resolution. 
Semiempirical two-phonon DOS spectra were used in Ref. 1^ to interpret the changes in 
Raman phonon line widths observed upon application of pressure and to assign structure 
in the two-phonon Raman spectra. We present here, for the sake of completeness, such 
calculations performed for the ab initio band structure of Fig. [H We display in Fig. [5]^a) 
the DOS that corresponds to the sums of two phonons (with total k=0), normalized to 36 
states per primitive cell. In order to identify the origin of the various bands involved, we 
also plot the integral of that DOS vs. frequency from uji + 002 equal to zero to its maximum: 
each plateau in the integrated DOS signals the end of a 2Ph band. These bands, and their 
weights, are labeled in the integrated curve. As an example we mention that the lowest band 
(from to ~ 220 cm~^) corresponds to two TA phonons, whereas the highest, rather narrow 
band corresponds to two LO phonons. We have plotted at the top of Fig. [5] the Raman 
spectra measured in this 2Ph region.- The various bands observed can, in this manner, be 
assigned to 2Ph combinations. The vertical (blue) lines try to establish correspondences 
between structure in the Raman spectra and in the 2Ph DOS. A more detailed assignment 

n 

to van Hove singularities (critical points) has been attempted in Ref|9|. 

Similar plots are shown in Fig. [6] for the calculated 2Ph difference DOS. The correspond- 
ing Raman spectra vanish at low temperatures because of the appropriate Bose-Einstein 
statistical factor (or, in other words, because no phonons are available for being annihi- 
lated). The spectrum on the top of Fig. [5](b) was obtained at 350 K, a temperature at 
which annihilation of TA phonons becomes possible. In this manner, the strong peaks at ~ 
220 cm~^ was identified as corresponding to the creation of a TO and the annihilation of a 
TA phonon. In spite of the additional structure in the 2Ph difference DOS shown in Fig. 
[HI no structure in the Raman spectra has been unambiguously identified as belonging to 
difference modes. We believe that the DOS of Fig. El combined with detailed measurements 
of the temperature dependence of the Raman spectra, should help in identifying additional 
2Ph difference structure. 
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FIG. 5: (color online) The lower vignette (c) displays the density of states which corresponds to 
the sum of two phonons with equal but opposite wavevector (2Ph DOS). The middle vignette 
(b) represents the integral of that 2Ph DOS (see text) whereas the upper vignette (a) shows the 
corresponding 2Ph Raman spectrum.— The 2Ph DOS is normalized to 36 states per primitive cell. 

V. SPECIFIC HEAT 



Several publications have already presented partial aspects of the specific heat of zb-ZnS. 
Low temperature data were early reported by Clusius and Harteck (19 K < T < 196 K) and 
by Martin (4 K < T < 20 K).-^"^ More recent data were presented by Birch but only for T < 
12 K.-^ Most of the existent high temperature data (above room temperature up to 1200 
K) were obtained for polycrystalline material of the type used in ir-windows (Cleartran® 
multispectral zinc sulfide).— Here we present experimental results for Cp{T) of four zb-ZnS 
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FIG. 6: Similar to Fig. [5]but for the DOS which corresponds to 2Ph differences. The normahzation 
is to 15 phonon difference states per primitive cell. 

samples grown by vapor phase transport (for details of the growth see Ref. 0) with different 
isotopic compositions: "^^^Zii^^^S (i.e. Zn and S with natural isotope composition), ^^Zn^^S, 
64Zn32S, 6^Zn34So.5 ^^So.s). The ^^Zn^^So.g sample should behave like ^^Zn^^S and thus, 

for simplicity, it was assumed to have this composition. The measurements covered the 
temperature range from 3 to 1 lOOK and so did also our corresponding ab initio calculations. 

The specific heat Cp between 2 and 280 K was measured with a PPMS systent^ as de- 
scribed before.-"^*^ Between room temperature and 1100 K the heat capacity was determined 
with a DSC 404 Fl Pegasus differential scanning calorimeter.— 

We show in Fig. [7] the Cp{T) data obtained for a sample with natural isotopic composition 
in our whole T range, together with our ABINIT-LDA calculations of C^{T). The difference 
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that appears above ~600 K, approximately linear in T, can be attributed to the thermal 
expansion and has been estimated with the standard expression:- 



C,{T) - C,{T) = al{T) B Kno; T (3) 

64, which in 



where a^iT) is the thermal volume-expansion coefficient taken from Ref. 
the temperature region under consideration increases slightly linear with the temperature. 
B is the isothermal bulk modulus and Vmoi the molar volume. Cp calculated from our 
using Eq. ([3]) is also shown in Fig. [71 It agrees rather well with the experimental results. 
For comparison, we have included the recently calculated ab initio by Hu et al. using 
the CASTEP code with the GGA approximation for the exchange and correlation potential 
at four temperatures.— These points agree rather well with our calculations. A small 
discrepancy appears at 300 K which could be due to the use of GGA vs. LDA: GGA leads 
to a larger lattice constant (Table [T]) and thus to lower phonon frequencies which correspond 
to higher values of 




r(K) 



FIG. 7: (color online) Measured (o and solid line) and calculated values of the specific heat of 
zb-ZnS in the 3-1200 K temperature range ((red) dotted line: ABINIT-LDA C^; (blue) dashed 
line: Cp calculated according to Eq. ([3]). The four points represented by (red) X correspond to 
the calculations of Hu et al.i^ The inset displays our high temperature data in an enlarged scale. 

The isotopic composition is not expected to significantly affect the high temperature 
results of Fig. [71 which must tend to the Petit and Dulong value of 49.9 J/molK.— The 
isotope effect should be seen best at the maximum in Cp/T^ which, for zb-ZnS is expected 
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to occur (and is in fact observed in the experimental data, cf. Fig. [8]) at around 21.5 K. 
This peak occurs in zb-type materials at a temperature that corresponds to ~l/6 of the 
frequency of the TA band in the DOS, expressed in temperature units, 130 K according to 
Fig. Ei 




FIG. 8: (color online) Temperature dependence of Cp/T^ in the T region from 3 to 50 K measured 
on a ZnS sample with natural isotope abundance. Low temperature data available in the literature 
are also plotted. The results of our ah initio ABINIT-LDA calculations are represented by the 
(red) solid line; extrapolation to T ^ K by the (red) dashed line. 

The measured effect of the substitution of ^^Zn for ^^Zn on the maximum of Cp/T^ 
vs. T is shown in Fig. [91 The Cp/T^ maximum is too shallow to be able to investigate the 
dependence of its temperature on isotopic mass, although it seems to shift up in temperature 
with decreasing mass, as expected. The value of Cp/T^ at its maximum, however, is shown 
to increase with increasing mass. This increase is nicely reproduced by our ABINIT- LDA 
calculations also shown in Fig. [9l 

Before proceeding with the isotopically modified samples we would like to recall that 
although most of the calculations presented here have been performed with the ABINIT- 
LDA code, we have also used the VASP code with either GGA or LDA exchange-correlation 
potentials. In order to compare the C„(T) results obtained with the three codes we display 
in Fig. [10] the corresponding computations for zb-ZnS with the natural isotopic abundances. 
The difference between the two LDA calculations ( ~3%) is not to be taken seriously, in 
view of the differences inherent to the VASP and ABINIT programs. 

More interesting is the higher value of the maximum in C„/T^ obtained from calculations 
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FIG. 9: (color online) Temperature dependence of Cp/T'^ in the T region from 3 to 50 K measured 
for two zb-ZnS samples with different isotopic compositions as indicated. Also given are the results 
of our ah initio ABINIT-LDA calculations. 
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FIG. 10: (color onUne) CJT'^ for nat2nnatg ^ calculated by ABINIT-LDA, and VASP-GGA and 
VASP-LDA, as indicated in the inset. 

using the GGA Hamiltonian. Since such calculations yield a larger lattice parameter than 
those based on LDA (see Table I), one can conjecture that the higher value of C^/T^ is 
related to the average lower phonon frequencies that result from the larger GGA lattice 
parameters. The measured values of Cy/T^ at the low temperature maximum fall between 
the LDA and the GGA calculations, as expected from the corresponding lattice parameters. 
The dotted values displayed in Fig. [10] below ~8 K have been drawn by free-hand: in this 
region the accuracy of our calculations does not allow reliable results. 

In Figure [TT] we display the logarithmic derivative of Cp/T'^ vs. the mass of the Zn isotope 
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^Zn evaluated for the experimental as well as the theoretical data of Fig. [91 Because of 
computational limitations, the calculated derivatives are meaningful only down to ~ 10 K. 
The dashed continuation of this curve to the value 1 has been performed by free hand. The 
value 1 is obtained using the expressions (derived from Debye's theory): 




FIG. 11: (color online) Logarithmic derivative of the specific heat versus Zn mass as measured and 
ab initio calculated by ABINIT-LDA for zb-ZnS. The dashed part of the calculated curve has been 
extrapolated free-hand so as to agree with the value of Eq.(|l]) at T = 0. 



dlnMzn 2 Mz„ + Ms ■ ^ ^ 

dHCjT^) ^ 3 Ms ^ Q ^ 

dlnMs 2 Mzn + Ms^ ' ^ ^ 

The maximum of the measurements of Fig. [H] takes place at 17 K. In most of our 
earlier work on this topic we have evaluated the ratio of the maximum frequency in the 
corresponding branch of the Ph-DOS to the temperature of the maximum in Cy/T^. In the 
case of Zn. the corresponding maximum is that of the TA phonons, which corresponds to 
130 K. The ratio (130 K/17 K) = 7.6 is rather close to that found for ZnO (7.3, Ref. U, 
a value which is also close to that calculated on the basis of a single frequency Einstein 
model.— 

Figure fT2] displays the logarithmic derivative oi C /T^ vs. the mass of the sulfur isotope, 
evaluated from experimental as well as theoretical data. Again, the theoretical curve has 
been extrapolated free-hand in the low temperature region so as to agree at T = with 
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Eq. ([5]). The calculations, as well as the experimental results, show two maxima which 
can be represented by two Einstein oscillators. The lowest maximum corresponds to TA 
phonons and its appearance reflects the existence of a considerable Zn component in the TA 
vibrations. The maximum at about 90 K can be assigned to the optical phonons (average 
frequency ~ 330cm~^ ~ 478 K). The ratio 478 K / 90 K = 5.3 is also close to that found 
for the corresponding maximum in ZnO (5.1)^). 
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FIG. 12: (color online) Logarithmic derivative of the specific heat versus S mass as measured and 
ab initio calculated by ABINIT-LDA for zb-ZnS. The dashed part of the calculated curve has been 
extrapolated free-hand so as to agree with the value of Eq.© at T = 0. 

In previous work on monatomic crystals we found a connection between the logarithmic 
derivative of C^/T^ versus T and the corresponding derivative vs. isotopic mass.- For binary 
materials, a similar connection was shown to hold provided one adds the two derivatives with 
respect to each of the isotope masses.— Figure [13] shows that this connection is also valid 
for ZnS. The experimental points in the flgure were obtained using the calculated mass 
derivatives and the calculated as well as the measured derivatives vs. T. 

VI. PRESSURE EFFECTS ON THE STRUCTURE AND LATTICE DYNAMICS 

The availability of reliable ab initio electronic structures enables us to tackle a few out- 
standing questions concerning the total energy and the lattice dynamics as a function of 
pressure. The flrst such question appeared in Sect. Ill in connection with the spin - orbit 
splitting at the top of the valence band in zb-ZnS. This splitting was calculated ab initio to 
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FIG. 13: (color online) Relationship between the derivative of C^,/T'^ vs. temperature and the 
corresponding derivatives with respect to the isotopic masses of Zn and S. The dotted lines were 
calculated from the heat capacities by varying the isotope mass of either Zn or S using the ABINIT- 
LDA code. The points are derived from the measured Cp/T'^ vs. T by taking numerical derivatives 
with respect to the isotope mass or with respect to temperature. 



be Ao = 0.062 eV, in reasonable agreement with the measured value (0.067 eV., Ref. |52| ). 
It was pointed out that this splitting was considerably smaller than that predicted on the 
basis of the 2p wave functions of S and the 3p of Zn. The reduction of s-o splitting was 
conjectured to be due to the admixture of ?)d core-like wave functions of Zn (this admixture 
is also responsible for the linear k-terms of Eq. ([1])). In order to test this conjecture we 
calculated the dependence of Aq on volume, which we express as: 

dlnAo 

For most tetrahedral semiconductors 7a is positive and lies around O.5.— This fact can 
be qualitatively understood as resulting from the wave function renormalization that takes 
place when the volume of the unit cell is changed. A uniform renormalization would lead 
to 7a =+1. The resistance of the core to the penetration of the wavefunction reduces this 
value to the typical values around +0.5 just mentioned. For ZnS, however, a change in 
volume is expected to alter the Srf-Zn admixture: a decrease in the volume should increase 
the 3d admixture at the top of the valence band. Because of the negative contribution of 
3d electrons to Aq this effect should result in a decrease in 7a, and possibly even a sign 
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reversal. Our ab initio VASP calculations yielded indeed 7a = -0.58. The value 7a = -0.12 
had already been obtained by Cerdeira et alM- from non-selfconsistent (not truly ab initio) 
Korringer-Kohn-Rostocker (KKR) calculations. 

The p-d admixture responsible for the anomalies just mentioned should not take place 
in crystals with inversion symmetry, e.g. in the rock salt phase of ZnS (rs-ZnS). We thus 
performed band structure calculations for rs-ZnS and obtained for p=0 (the rs-ZnS phase 
is not stable at this pressure but this should not affect our argument) 7a = -1-0.60 thus 
confirming our conjecture concerning the nature of Aq in both phases of ZnS. As already 
mentioned in Sect. IIIIl the value of Aq calculated for the rs phase is 0.112 eV, much larger 
than that of the zb phase. 

In order to proceed with our studies of the effects of pressure on the electronic and 
vibrational properties of ZnS we have investigated the region of stability of the zb and also 
the wurtzite (wz) phase concerning the transition to the rs phase.— 

A number of calculations of the pressure at which the zb to rs transition and also the wz 
to rs transition take place have already appeared in the literature. Most of them use the 
LDA Hamiltonian^i^'^'^'^, one of them uses the GGA Hamiltonian.— The extant LDA 
calculations indicate transition pressures from zb-ZnS to rs-ZnS between 13 and 19.5 GPa, 
whereas the GGA calculation gives 17.2 GPA.-- Experimental values lie between 14.7 and 
18.1 GPa.— Because of the broad range of transition pressures which have been reported, 
we decided to perform calculations using both, LDA and GGA Hamiltonians. The obtained 
enthalpy differences between rs, wz and the zb phase are shown in Fig. [TUfor the VASP-LDA 
calculation (including s-o interaction). In this figure the zero temperature enthalpy of the 
wz modification lies about 0.008 eV per primitive cell above that of zb (it decreases slightly 
towards zero pressure) The zb to rs transition is predicted to occur at 15.7 GPa, a pressure 
that falls in the middle of the thus far reported range (13 to 19.5 GPa). Our corresponding 
VASP-GGA calculation yielded a wz phase also 0.008 eV per primitive cell above zb and 
a transition pressure of 16.7 GPa. One can also here conjecture that the higher transition 
pressure of the GGA calculation is related to the larger lattice parameter. 

In order to investigate additional effects of the inversion symmetry present in the rs 
structure we calculated the electronic band structure of rs-ZnS at a pressure p = 16.0 GPa, 
at which this phase is stable (see Fig. [T^ . The results are displayed in Fig. [151 The s-o 
splitting of the valence bands at the Fpoint is found to be 0.12 eV, considerably higher than 



23 



0.6 



0.4 



> 
a. 

CD 



0.2 



c 
LU 



-0.2 



-0.4 



■ 


ZnS (LDA SO) 


wurtzite 


\ P 


j= 15.7 GPa 


zb 












. 1 , 


1 1 


. 1 



5 10 15 20 
Pressure (GPa) 



25 



FIG. 14: Enthalpy difference between the rs and wz phases of ZnS and that of the zb phase (the 
latter is assumed to have zero enthalpy) at T= K (neglecting zero-point vibrations), as calculated 
with the VASP-LDA code. 

that of the zb modification, as already discussed. The most conspicuous difference between 
the band structure of Fig. [15] and that of Fig. [1] is the lack of a gap: The minimum of 
the conduction band is found at the X-point, at an energy of -0.2 eV. The maximum of the 
valence bands is along the F-K direction, at an energy of 0.1 eV. Thus there seems to be 
an overlap of about 0.3 eV between valence and conduction bands. At this point, however, 
we must keep in mind the "gap problem" already mentioned in Sect. IIIII For the case 



of the zb-ZnS modification, the calculated va^ 
apart by about 2.05 eV (Our Refs. 



47 and 



ence and conduction bands had to be pulled 



571 ) in order to correct this problem. Using the 



same correction for the band structure of Fig. [T5] as a so-called " scissors operator" , a gap 
of ~1.75 eV opens and rs-ZnS becomes an indirect gap semiconductor, with the valence 
band maximum along F-K (S direction) and the conduction band minimum at X. A band 
structure calculation for rs-ZnS has already appeared in the literature.— It was performed 
with the linear muffin tin orbital method (LMTO) also using an LDA Hamiltonian and was 
thus also affected by the gap problem. The authors of Ref. l72 performed transmission 
measurements for a rs-ZnS crystal obtained under pressure. They found an indirect gap of 
~1.9 eV, in good agreement with the one we obtained after applying the "scissors operator" 
as mentioned above 
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FIG. 15: Band structure of rs-ZnS calculated for p=16 GPa using the ABINIT-LDA code (s-o 
coupling included). Notice the overlap between conduction and valence bands which, as explained 
in the text, is spurious and disappears once the appropriate "scissors operator" is applied. 

It is easy to explain in a qualitative manner the remarkable difference between the band 
structures of Figs. [15] and [1] as due to the effect of inversion symmetry (which is lacking for 
zb). Inversion symmetry prevents the hybridization of the p-like top of the valence bands 
with the 3d core states at the F point of zinc since both sets of states have opposite parity. 
When moving away from the F point, either towards L or K, parity is no longer a good 
quantum number, the p and d states mix and the p-like top of the valence band is pushed 
up, as shown in Fig. [151 This suffices to make the gap indirect. 

The effect just described is similar to that encountered in rs-structure AgCl and AgBr 
(Ref. QS) which also have a valence band maximum along S and an indirect gap around 
3 eV. The minimum of the conduction band of these materials, however, is at the Fpoint 
at p=0. Application of pressure (our calculations were per 



brmed for p=16 GPa) would 
probably raise the F minimum and lower that at X (Ref. |75|) bringing the band structure 
of AgCl and AgBr even more in line with that of Fig. [151 

We conclude this section by showing calculations of the pressure dependence of the 
phonons of natural ZnS calculated with the ab initio ABINIT-LDA code at several high 
symmetry points (F, X and L). The figure displays the calculated points (full circles ) and 
also experimental points obtained by Weinstein.— The solid lines are quadratic fits to the 
calculated results. Of particular interest is the supralinear decrease of the TA frequencies 
with pressure found at the edge of the Brillouin zone (X and L points). This decrease 
seems to be qualitatively related to the rs phase transition^ The agreement between the 
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experimental and calculated results of Fig. [16] is remarkably good. 




FIG. 16: (solid lines) Pressure dependence of the phonons of natural ZnS as calculted with the 
ab initio ABINIT-LDA code at the indicated high symmetry points in the Brillouin zone. Experi- 
mental data taken from Ref. |77|. are represented by x 



From the coefficient of the linear terms of the fits in Fig. [16] we have calculated the 
corresponding mode Griineisen parameters using the calculated bulk modulus Bq = 85.86. 
They are listed in Table [TT] together with extant experimental data. Again, the agreement 
is good. Notice that at the 7 point 7x0 is larger than 7lO) a fact common to most zb-zype 
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semiconductors.— These Griineisen parameters can be used to calculate the dependence of 
the transverse effective charge on volume, and the corresponding Griineisen parameter 7e, 
using the expression obtained from Eq. ([2]): 



dlnV 2dlnV uj'^q — uj\q 

The logarithmic derivative of for ZnS has been reported by Vedam and Schmidt to 
be +0.36.— Using this value and those from Table [TTl into Eq. ([7j) we find 7e = -0.84, a 
result rather similar to that found for other zb-type semiconductors (exception: SiC).— 
The negative value of 7e implies an increase in ionicity with increasing volume. Using the 
experimental values of the 7 and u in Table HTl we obtain 7e = -0.44. 



TABLE II: Mode Griineisen parameters at the F, L and X point as obtained from the pressure 
dependence of the phonon energies obtained by ABINIT-LDA calculations compared with avail- 
able experimental data. If we use the experimental value Bq= 76.9 instead of the ABINIT-LDA 
calculated one {Bq= 85.86) the agreement between the calculated and the experimental values 
improves) . 
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VII. CONCLUSIONS 



The original motivation of this work was to perform accurate measurements of the specific 
heat of zb-ZnS over a broad temperature range, using samples of different isotopic composi- 
tions and to compare these with ab initio calculations based on phonon dispersion relations 
obtained from ab initio band structures. Because of the availability of several band structure 
codes, we also decided to compare the results obtained using them with either LDA or GGA 
exchange-correlation Hamiltonians. The agreement between experimental and theoretical 
results has been found to be quite satisfactory, the former lying between the rather close 
LDA and GGA calculations. Because of the extensive computational tools we developed, 
we decided to use them to tackle a few extant problems involving the electronic structure 
and the lattice dynamics. We considered the effect of spin-orbit interaction on the lattice 
parameters and the lattice dynamics and found it to be negligible (contrary to the results 
obtained earlier for materials with heavier atoms, such as PbTe or Bi). We conjectured that 
the anomalously small s-o splitting at the top of the valence bands of zb-ZnS was due to 
admixture of Sd-Zn core levels. In order to prove this conjecture, we performed calculations 
for ZnS in the centrosymmetric rock salt structure. In this structure no 3d-Zn admixture is 
possible at the top of the valence bands and the s-o splitting has the larger expected value. 
We also calculated the pressure dependence of the s-o splitting of zb-ZnS and found it to 
anomalously decrease with pressure, a fact that we have attributed to the pressure depen- 
dent contribution of the 3d-Zn electrons. Another effect of the 3d admixture is reflected in 
the valence bands of rs-ZnS which, like those of similar rs-type materials (e.g. AgCl) have 
a maximum away from the center of the Brillouin zone (our Ref. |l9|. We have confirmed 
this by calculating the band structure of this phased which erroneously turned out to be 
semimetallic because of the "gap problem". By applying an ad hoc "scissors operator" we 
found that the indirect gap should be about 1.75 eV, in reasonable agreement with that 
found experimentally (1.9 eV, Ref. l72l ). The zb phase of ZnS has linear splittings of the va- 



lence bands around the F point of the BZ. We have calculated these splittings which should 
be of interest in the field of spintronics. 

We have also investigated pressure effects on the crystal structure of ZnS. The zb-rs phase 
transition was calculated to occur at 15.7 GPa when using the LDA Hamiltonian, whereas 
the use of the GGA led to a transition pressure about 1 GPa higher, a fact which is probably 
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related to the different lattice parameters found at zero pressure when using either of the 
two approximations. We have performed several studies of the lattice dynamics of zb-ZnS. 
The ab initio calculations agree rather well with the experimental data obtained by INS and 
Raman spectroscopy. We have also presented calculation of the phonon density of states 
(including its projection on each of the constituent atoms) and the optical two-phonon DOS. 
The latter has been compared with second order Raman spectra. 

Finally, we have calculated the pressure dependence of the phonon frequencies at the F, X 
and L high symmetry points and the corresponding Griineisen parameters. The Griineisen 
parameter of the transverse (Born) effective charge has also been evaluated. 
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